C----------------- Ux ------- UU(I,J)=UU(I,J)+U(I,J)*UX*H ------------------
DO 45 J=1,N
DO 43 I=2,N
B(I)=4.
A(I)=1.
C(I)=1.
S(I)=3*(U(I+1,J)-U(I-1,J))
43 CONTINUE
CALL SOLVER(A,B,C,S,G,H1,2,N)
DO 44 I=2,N
UXH1=S(I)
44 UU(I,J)=U(I,J)*UXH1
UXXHHL(J)=S(2)-0*H
UXXHHR(J)=0*H-S(N)
45 CONTINUE
C----------------- Vy -------- VV(I,J)=VV(I,J)+V(I,J)*VY*H -----------------
DO 355 I=1,N
DO 353 J=2,N
B(J)=4.
A(J)=1.
C(J)=1.
S(J)=3*(V(I,J+1)-V(I,J-1))
353 CONTINUE
CALL SOLVER(A,B,C,S,G,H1,2,N)
DO 354 J=2,N
VYH=S(J)
354 VV(I,J)=V(I,J)*VYH
VYYHHD(I)=S(2)-0*H
VYYHHU(I)=0*H-S(N)
355 CONTINUE
C------------ Uxx -------- UU(I,J)=UU(I,J)-UXX*H/RE ------- UA=UXXHH ----------
DO 65 J=1,N
DO 63 I=3,NM1
A(I)=1.
B(I)=10.
C(I)=1.
63 S(I)=(U(I-1,J)-2.*U(I,J)+U(I+1,J))*12.
IF(NORDER.LE.2) THEN
B(2)=1
C(2)=0
S(2)=U(1,J)-2*U(2,J)+U(3,J)
A(N)=0
B(N)=1
S(N)=U(NM1,J)-2*U(N,J)+U(N1,J)
ELSE
B(2)=16
C(2)=2
S(2)=U(1,J)*27-48*U(2,J)+21*U(3,J)
A(N)=2
B(N)=16
S(N)=21*U(NM1,J)-48*U(N,J)+U(N1,J)*27
END IF
CALL SOLVER(A,B,C,S,G,H1,2,N)
DO 64 I=2,N
UXXHH=S(I)
UA(I,J)=UXXHH
64 UU(I,J)=UU(I,J)-UXXHH/REH
65 CONTINUE
C---------------- Uyy ------- UU(I,J)=UU(I,J)-UYY*H/RE ----- UA ---------------
DO 365 I=2,N
DO 363 J=2,NM1
A(J)=1.
B(J)=10.
C(J)=1.
363 S(J)=(U(I,J-1)-2.*U(I,J)+U(I,J+1))*12.
B(1)=5-0.125
C(1)=1+0.25
S(1)=8*(U(I,2)-3*U(I,1)+U0*2)
A(1)=-0.125 ! COEFF OF Uyy(I,3)
C(N)=-0.125 ! COEFF OF Uyy(I,N-2)
A(N)=1+0.25
B(N)=5-0.125
S(N)=8*(U(I,NM1)-3.*U(I,N)+2*U1)
CALL SOLVER1(A,B,C,S,G,H1,1,N)
DO 364 J=1,N
UYYHH=S(J)
UA(I,J)=UA(I,J)+UYYHH
364 UU(I,J)=UU(I,J)-UYYHH/REH
DO 367 J=N,2,-1
367 UA(I,J)=(UA(I,J-1)+UA(I,J))/2
IF(I.EQ.2) THEN
DO J=1,N
UYYHHL(J)=(0*HH+S(J))/2.
ENDDO
ELSEIF(I.EQ.N) THEN
DO J=1,N
UYYHHR(J)=(0*HH+S(J))/2.
ENDDO
ENDIF
365 CONTINUE
C----------------- Vx ------- VV(I,J)=VV(I,J)+UA*VX*H ------------------
DO 55 J=2,N
IF(NORDER.LE.2) THEN
DO 3661 I=1,N
3661 UA(I,J)=0.25*(U(I,J)+U(I+1,J)+U(I,J-1)+U(I+1,J-1))
ELSE
UA(1,J)=0.25*(U(1,J)+U(2,J)+U(1,J-1)+U(2,J-1))
* -0.5*(UXXHHL(J-1)+UXXHHL(J)
* +UYYHHL(J-1)+UYYHHL(J))/8.
DO 366 I=2,NM1
366 UA(I,J)=0.25*(U(I,J)+U(I+1,J)+U(I,J-1)+U(I+1,J-1))
& -0.5*(UA(I,J)+UA(I+1,J))/8
UA(N,J)=0.25*(U(N,J)+U(N1,J)+U(N,J-1)+U(N1,J-1))
* -0.5*(UXXHHR(J-1)+UXXHHR(J)
* +UYYHHR(J-1)+UYYHHR(J))/8.
END IF
DO 53 I=2,NM1
B(I)=4.
A(I)=1.
C(I)=1.
S(I)=3*(V(I+1,J)-V(I-1,J))
53 CONTINUE
B(1)=9
C(1)=3
S(1)=8*(V(2,J)-V0)
A(N)=3
B(N)=9
S(N)=8*(V1-V(NM1,J))
CALL SOLVER(A,B,C,S,G,H1,1,N)
DO 54 I=1,N
VXH=S(I)
54 VV(I,J)=VV(I,J)+UA(I,J)*VXH
55 CONTINUE
C---------------- Vxx ------- VV(I,J)=VV(I,J)-VXX*H/RE -------------------
DO 75 J=2,N
DO 73 I=2,NM1
A(I)=1.
B(I)=10.
C(I)=1.
73 S(I)=(V(I-1,J)-2.*V(I,J)+V(I+1,J))*12.
B(1)=5-0.125
C(1)=1+0.25
S(1)=8*(V(2,J)-3*V(1,J)+V0*2)
A(1)=-0.125 ! COEFF OF Vxx(3,J)
C(N)=-0.125 ! COEFF OF Vxx(N-2,J)
A(N)=1+0.25
B(N)=5-0.125
S(N)=8*(V(NM1,J)-3.*V(N,J)+V1*2)
CALL SOLVER1(A,B,C,S,G,H1,1,N)
DO 74 I=1,N
VXXHH=S(I)
VA(I,J)=VXXHH
74 VV(I,J)=VV(I,J)-VXXHH/REH
IF(J.EQ.2) THEN
DO I=1,N
VXXHHD(I)=(0*HH+S(I))/2.
ENDDO
ELSEIF(J.EQ.N) THEN
DO I=1,N
VXXHHU(I)=(0*HH+S(I))/2.
ENDDO
ENDIF
75 CONTINUE
C---------------- Vyy --------- VV(I,J)=VV(I,J)-VYY*H/RE -----------------
DO 375 I=1,N
DO 373 J=3,NM1
A(J)=1.
B(J)=10.
C(J)=1.
373 S(J)=(V(I,J-1)-2.*V(I,J)+V(I,J+1))*12.
IF(NORDER.LE.2) THEN
B(2)=1
C(2)=0
S(2)=V(I,1)-2*V(I,2)+V(I,3)
A(N)=0
B(N)=1
S(N)=V(I,NM1)-2*V(I,N)+V(I,N1)
ELSE
B(2)=16
C(2)=2
S(2)=V(I,1)*27-48*V(I,2)+21*V(I,3)
A(N)=2
B(N)=16
S(N)=21*V(I,NM1)-48*V(I,N)+V(I,N1)*27
END IF
CALL SOLVER(A,B,C,S,G,H1,2,N)
DO 374 J=2,N
VYYHH=S(J)
VA(I,J)=VA(I,J)+VYYHH
374 VV(I,J)=VV(I,J)-VYYHH/REH
375 CONTINUE
DO 376 I=N,2,-1
DO 376 J=2,N
376 VA(I,J)=(VA(I-1,J)+VA(I,J))/2
C----------------- Uy ------ UU(I,J)=UU(I,J)+VA*UY*H ---------------------
DO 345 I=2,N
IF(NORDER.LE.2) THEN
DO 3781 J=1,N
3781 VA(I,J)=0.25*(V(I,J)+V(I-1,J)+V(I,J+1)+V(I-1,J+1))
ELSE
VA(I,1)=0.25*(V(I,1)+V(I-1,1)+V(I,2)+V(I-1,2))
& -0.5*(VXXHHD(I-1)+VXXHHD(I)
* +VYYHHD(I-1)+VYYHHD(I))/8
DO 378 J=2,NM1
378 VA(I,J)=0.25*(V(I,J)+V(I-1,J)+V(I,J+1)+V(I-1,J+1))
& -0.5*(VA(I,J)+VA(I,J+1))/8
VA(I,N)=0.25*(V(I,N)+V(I-1,N)+V(I,N1)+V(I-1,N1))
& -0.5*(VXXHHU(I-1)+VXXHHU(I)
* +VYYHHU(I-1)+VYYHHU(I))/8
END IF
DO 343 J=2,NM1
B(J)=4.
A(J)=1.
C(J)=1.
S(J)=3*(U(I,J+1)-U(I,J-1))
343 CONTINUE
B(1)=9
C(1)=3
S(1)=8*(U(I,2)-U0)
A(N)=3
B(N)=9
S(N)=8*(U1-U(I,NM1))
CALL SOLVER(A,B,C,S,G,H1,1,N)
DO 344 J=1,N
UYH=S(J)
344 UU(I,J)=UU(I,J)+VA(I,J)*UYH
345 CONTINUE
IF (KSTEP.EQ.1) THEN
IF (KU.EQ.1) THEN
DO 350 I=2,N
DO 350 J=1,N
UUN(I,J)=UU(I,J)
350 VVN(J,I)=VV(J,I)
ELSE ! KU=2
DO 352 I=2,N
DO 352 J=1,N
UUN(I,J)=UUN(I,J)+2*UU(I,J)
352 VVN(J,I)=VVN(J,I)+2*VV(J,I)
END IF
ELSE IF (KSTEP.EQ.2) THEN
DO 360 I=2,N
DO 360 J=1,N
UUN(I,J)=UUN(I,J)+2*UU(I,J)
360 VVN(J,I)=VVN(J,I)+2*VV(J,I)
ELSE IF (KSTEP.EQ.3) THEN
DO 362 I=2,N
DO 362 J=1,N
UU(I,J)=(UUN(I,J)+UU(I,J))/6
362 VV(J,I)=(VVN(J,I)+VV(J,I))/6
ENDIF
C---------------- Px ------ U(I,J)=UST(I,J)-PX*DT -------------------
DO 85 J=1,N
DO 83 I=3,NM1
A(I)=1
B(I)=22
C(I)=1
S(I)=24*(P(I,J)-P(I-1,J))
83 CONTINUE
B(2)=25
C(2)=-2
A(2)=1
S(2)=24*(P(2,J)-P(1,J))
C A(2) IS COEFF OF Px(4,J), C(N) IS COEFF OF Px(N-2,J)
C(N)=1
B(N)=25
A(N)=-2
S(N)=24*(P(N,J)-P(NM1,J))
CALL SOLVER1(A,B,C,S,G,H1,2,N)
DO 81 I=2,N
PXH=S(I)
81 U(I,J)=UST(I,J)-PXH*DTDH2
85 CONTINUE
C---------------- Py ---------- V(I,J)=VST(I,J)-PY*DT ------------------
DO 385 I=1,N
DO 383 J=3,NM1
A(J)=1
B(J)=22
C(J)=1
S(J)=24*(P(I,J)-P(I,J-1))
383 CONTINUE
B(2)=25
C(2)=-2
A(2)=1
S(2)=24*(P(I,2)-P(I,1))
C A(2) IS COEFF OF Px(I,4), C(N) IS COEFF OF Px(I,N-2)
C(N)=1
B(N)=25
A(N)=-2
S(N)=24*(P(I,N)-P(I,NM1))
CALL SOLVER1(A,B,C,S,G,H1,2,N)
DO 384 J=2,N
PYH=S(J)
384 V(I,J)=VST(I,J)-PYH*DTDH2
385 CONTINUE
C----------------- Central Ux*H ------ UXH ---------------------
DO 845 J=1,N
DO 843 I=2,NM1
A(I)=1
B(I)=22.
C(I)=1
S(I)=24*(U(I+1,J)-U(I,J))
843 CONTINUE
B(1)=15
C(1)=1
S(1)=18*(U(2,J)-U(1,J))
A(N)=1
B(N)=15
S(N)=18*(U(N1,J)-U(N,J))
CALL SOLVER(A,B,C,S,G,H1,1,N)
DO 844 I=1,N
844 UXH(I,J)=S(I)
845 CONTINUE
C----------------- Central Vy*H ----- DIVH -------------------
DO 855 I=1,N
DO 853 J=2,NM1
A(J)=1
B(J)=22
C(J)=1
S(J)=24*(V(I,J+1)-V(I,J))
853 CONTINUE
B(1)=15
C(1)=1
S(1)=18*(V(I,2)-V(I,1))
A(N)=1
B(N)=15
S(N)=18*(V(I,N1)-V(I,N))
CALL SOLVER(A,B,C,S,G,H1,1,N)
DO 854 J=1,N
VYH=S(J)
854 DIVH(I,J)=UXH(I,J)+VYH
855 CONTINUE
C------------- SUBROUTINE SOLVER -------------------
SUBROUTINE SOLVER(A,B,C,S,G,H,N0,N)
to get solution of
A(i)*X(i-1)+B(i)*X(i)+C(i)*X(i+1)=S(i), (i=N0,...,N)
X(N0-1)=X(N+1)=0
then set S(i)=X(i), (i=N0,...,N)
C------------- SUBROUTINE SOLVER1 -------------------
SUBROUTINE SOLVER1(A,B,C,S,G,H,N0,N)
the difference between "SUBROUTINE SOLVER" and this SUBROUTINE is:
the equations for i=N0 and i=N are
B(i)*X(i)+C(i)*X(i+1)+A(i)*X(i+2)=S(i), (i=N0)
C(i)*X(i-2)+A(i)*X(i-1)+B(i)*X(i)=S(i), (i=N)
(if set A(N0)=C(N)=0, then these two SUBROUTINEs are same).
C----------------------- SUBROUTINE MG -------------------------
SUBROUTINE MG(P,DP,P1,PA,NA,R,DR,R1,RA,A,B,C,S,G,H1,SS
& ,N,N1,AKPMAX,ERPL,EPSP,KQMG,EPSPA,KQW,PP)
MultiGrid method